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ABSTRACT 

Several integration schemes exits to solve the equations of motion of the iV-body 
problem. The Lie-integration method is based on the idea to solve ordinary differen- 
tial equations with Lie-series. In the 1980s this method was applied for the iV-body 
problem by giving the recurrence formula for the calculation of the Lie-terms. The 
aim of this works is to present the recurrence formulae for the linearized equations 
of motion of iV-body systems. We prove a lemma which greatly simplifies the deriva- 
tion of the recurrence formulae for the linearized equations if the recurrence formulae 
for the equations of motions are known. The Lie-integrator is compared with other 
well-known methods. The optimal step size and order of the Lie-integrator are calcu- 
lated. It is shown that a fine-tuned Lie-integrator can be 30%-40% faster than other 
integration methods. 
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1 INTRODUCTION 

The classical problems of celestial mechanics are described 
by a system of ordinary differential equations (ODEs). The 
investigation of the motions in the Solar System, exoplan- 
etary systems, satellites around the Earth or other celes- 
tial objects are based on the solutions of such ODEs. How- 
ever, several modern analysis, including many chaos detec- 
tion methods require to solve the linearized equations of the 
problem. 

The integration m ethod based on the Lie-series 
l|Grobner fc Knapdll967f) is widely used in celestial m echan- 
ics to solve ODEs (see Hanslmeier fc Dvorakl l| 19841 ). here- 
after H&D and articles referring to it). The basis of this 
method is to generate the coefficients of the Taylor expan- 
sion of the solution by using recurrence relations. The prin- 
cipal application, i.e. the integration of the iV-body problem 
is described in details in H&D. 



1.1 Lie-integration 

Here we summarize the key points of this method of numer- 
ical integration, using almos t identical notations as used by 
lHanslmeier fe Dvorakl (|l984f ). 

Let us write the differential equation to be solved as 

*< = /<(x), (1) 
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where x = (xi, . . . , xn) is an R — + R^ and f = (/i, . . . , /jv) 
is an — > R" continuous function and N is the dimension 
of the vector x and the vector space where f maps from and 
maps to. Let us introduce the differential operator 



Di := 



d 



dxi 

and the derivation 



d 

dxi ' 



(2) 



(3) 



which is known as the Lie-derivation or Lie-operator. Lo is a 
linear differential operator and one can apply Leibniz's rule, 



L (ab) = aL (b) + bL (a) 



(1) 



where a and b are R — > R differentiable functions. It can 
easily be proven that the solution of equation (TTJ at a given 
instance t + At is formally 



x(t + At) = exp (At ■ L ) x(t) 
where 

At k 



exp(At-Lo) = Y, 



A:! 



-L 



fc! 



(5) 



(6) 



The method of Lie-integration is finite approximation of the 
sum in the right-hand side of equation ((6j), up to the order 
of M, namely 
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x(i + At) 



\ -» At k 



^) = J2^?( L "Mt))- (7) 



The proof of equation ([5]) and other related proper t ies of the 
Lie-derivation can be f ound in iGrobner fc Knappl (|l967T ) or 
iHanslmeier fc Dvorakl (|l984T h 

In spite of the fact that the Lie-derivatives can ana- 
lytically be calculated up to arbitrary order, the formulae 
yielded by these expansions are highly complicated even if 
all kinds of new variables are introduced (see e.g. equations 
(19d) or (19e) in H&D at page 204). A definitely more effi- 
cient way to evaluate the Lie-derivatives is to find a set of 
recurrence relations. These relations allow us to express the 
(n + l)th Lie-derivative, e.g. L n+1 x as the function of the 
derivatives with lower order, namely L J x where ^ j ^ n. 
The initialization of such a recurrence relation is evident, 
because L°x = x. We note that in several applications, well- 
chosen auxiliary variables have to be introduced to gain a 
compact set of recurrence relations which can efficiently be 
evaluated. 



1.2 The importance of linearized equations 

Wide range of problems related to celestial mechanics re- 
quire to solve simultaneously the linearized form of the orig- 
inal equations too. The numerous experiments conducted in 
the last decades show that chaotic behaviour is typical and 
already occurs in simple but nonlinear systems. This find- 
ing throws completely new light upon these systems and the 
study of chaotic behaviour became of high concern. A major 
part of the frontline research focuses on the structure of the 
phase space, therefore the problem to separate ordered and 
chaotic motion in systems, which posses only a few degrees 
of freedom and are described by ODEs, has become a fun- 
damental task in a wide area of modern research. The phase 
space of these nonlinear systems can not be described by 
the known mathematical tools. To map the phase space and 
study the chaotic behaviour of a given system fast and reli- 
able numerical tools are needed. These tools are extremely 
useful in those cases when the inspected dynamical system 
has more than two degrees of freedom and accordingly its 
phase space cannot be explored in a direct way or the clas- 
sical method of surface of section (SoS) can not be applied 
which is widely used in the case of conservative systems with 
two degrees and fre edom. The basi c idea of the method of 
SoS was inv ented bvlPoincare (I1899T ) and its application was 
renewed bv lHenon fc Heilesl 1 1964 ). 

The mathematical foundation of the theory of Lya- 
punov Characteristic Exponents (LCEs) is approximately of 
the same age as the SoS and arose progressively in the lit- 
eratur e. The use of such expone nts dates b ack t o iLvapunovl 
(|l907h . but was fir st applied bylOseleded |l968) to charac- 
terize trajectories. iHenon fc Heilesl jj 19641 ) found that in an 
integrable region of the phase space of a dynamical system 
nearby orbits diverge linearly whereas in a chaotic region 
they diverge exponentially. The LCEs express these facts in 
a precise form and many papers were devoted to the appli- 
cation of LCEs in several nonlinear problems. 

Unfortunately both methods have a serious drawback. 
To compute the LCEs the equations have to integrate for in- 
finity, which is numerically impossible. The method of SoS 
becomes hard to handle and greatly deceiving for systems 



with more than two degrees of freedom. To overcome these 
problems was the main motivation in the 1980s that initiated 
the research to develop new numerical methods to character- 
ize the stochasticity of the trajectories in the phase space in 
short time-span and in arbitrary dimension. The developed 
methods can be classified in two groups: one group con- 
sists of the methods which are based on th e analys i s of t he 
orbits, (e.g. SoS or frequency analysis, see iLaskarl |l990)), 
the other one is based on the time evolution of the tangent 
vector, i.e. the solution of the linearized equations of mo- 
tion (e.g. LCE). There are complete software packages de- 
signed to analyse systems of celestial mechanics, both for for 
gener al integration of motion (e.g. Mercury6, see IChambersI 
Il999f ) and for solvin g linearized equations and calculating 
LCEs (ORBIT9, see llVlilani fc Nobilil IT98SI ) . We also have 
to mention that there are several improved chaos detec- 
tion methods which are based on the solution of the lin- 
earized equations. Instead of a complete review, we only 
mention two of t hem: the method of Fast Lyapunov Indi- 
cators (FLIs, see iFroeschle et al1ll997f ) and the method of 
Mean Exponential Growth of Nearby Orbits (M EGNO, see 
ICincotta fc Sim6ll2000l : iGodziewski et al.ll200lf ). 

The aim of this paper is to present a lemma which ad- 
vances the derivation of the same kind of recurrence relations 
for the linearized equations. We present these relations for 
certain classical dynamic systems: for the general A-body 
problem and for the A-body problem in the reference frame 
of one of the bodies. In the last section we compare the 
efficiency of this method with well-known other ones. 



2 LINEARIZED EQUATIONS 

The chaos indicators mentioned in the previous section can 
be obtained if the linearized form of the equations of motion 
is solved. The solution of the linearized equations is an £ = 
■ ■ ■ , £.n) : R — > function having the same dimension 
as the equations of motion has. The linearized equations of 
equation (QJ (any ODE can be written in this form) can be 
written as 



OX m 



(8) 



where the variables £j = £i(t) : R — > R are the so-called 
linearized variables and x = x(t) is the solution of equa- 
tion |l}. Equation ([8]) is linear in £, therefore if Q (t) and 
£j (t) are two independent solutions, then a^ 1 ' (t) +f3(,f^ (t) 
is also a solution. Using the Einstein summation convention 
equation (|8} can be written in a more compact form: 

£,i ^mD m fi. (9) 

2.1 Lie-derivatives of the linearized equations 

Introducing the differential operator 

d > ■= w (10) 

the coupled system of equations (both the original and the 
linearized) is 

Xi = ft, (11) 

ii = £mD m fi, (12) 
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and the Lie-operator of equations (|11[) - (|12[ ) i s 

L = L + L t = f l D, + i m D m f t d l . (13) 

Lemma Using the same notations as above the Lie-deriva- 
tives of £fc can be written as 



L n £k — £,mDmL n Xk — (,mD m LQXk- 



(14) 



Proof Obviously, equation (|14[) is true for n — 0: 
D m L°Xk = D m Xk = 5 m k, (15) 
hence 

£mDmL°Xk = £,m5mk = £fc ■ (16) 

Let us suppose that it is true for all ^ j ^ n and calculate 
the (n + l)th Lie-derivative of 

L n+1 6 = L (UDmL n x k ) = 

= (/*A + ijDjfidi) (^ m D m L n x k ) = 
= fiDi£mDmL n x k + 

+£i(Djfi)[5imD m L"x k + 

+£™D m diL"x k }. (17) 

Here the term £,mDmdiL n x k equals to zero, because x k and 
L n Xk for all n ^ do not depend on £. Therefore, 

= f 1 D l £ m D m L n x k + £j(Djfi)DiL n Xk = 

= £ m fiDmDiL n x k + t; m (Dmfi)(DiL n x k ) = 

= £m (fiD m + D m fi) (DiL"x k ) = 

= imD m {fiD l )(L n x k ) = 

— ^mDmL^L x k ) = ^ m D m L x k 

= UDmL^Xk. (18) 

We have applied Young's theorem, namely 

D m Di = DiD m , (19) 

and Leibniz rule, 

D m {fiDi)X = D m fi(DiX) = 

= fi{D m DiX) + {D m fi){DiX), (20) 

where X can be an arbitrary function of x, in equation (|18p 
X = L n x k . Therefore equation (|18fl is the same relation 
for n + 1, as equation (|14[1 for n. Continuing the scheme 
described above, equation (I14|l can be proven for all positive 
integer values of n. ■ 

2.2 An example: applying to the Henon-Heiles 

system 

Demonstrating the power of the lemma proven in the pre- 
vious subsection we derive the recurrence relations for the 
equation of the Henon-Heiles dynamical system and its lin- 
earized form. The Henon-Heiles system is one of the simplest 
Hamiltonian systems which s hows chaotic behaviou r under 
certain initial conditions (see lHenon fc Heileslll964h . 

The equations of motion are derived from the Hamilto- 
nian function 

H(x,y;v,w) = i (^x 2 + y 2 + 2x 2 y - ^y 3 ^j + 

+ 1 (y 2 +W 2 ) , 



where x — v and y = w. The equations of motion are 

x = v, (22) 

y = w, (23) 

v = —x — 2xy, (24) 

w = -y-x 2 + y 2 , (25) 

and the Lie-operator of this system of equations is 

Lo = vd x + wd y + (-x — 2xy)d v + (—y - x 2 + y 2 )d w , (26) 

according to equation J3J. It can easily be shown that the 
recurrence relations of the equations (|22|l - (I25p are the fol- 
lowing, 

Lq +1 x = LqV, (27) 
L n +l y = L%w, (28) 



Ln + 1 j n n\ ^ 1^1 rk T n — k 

o v = -L x~2 2_^ \L xL y, 



L ° y ~ ( fe ( L o xL o 



oyL y 



(29) 



(30) 



Let us denote the linearized variables related to x, y, v and 
w by £, rj, 4> and p, respectively. According to equation (|14|l 
the Lie-derivatives of these variables are 



L n £ = (, m D m L n x, 

L n rj = £mD m L n y, 

L n 4> — (, m D m L n v, 

L n p = £ m D m L n w, 



(31) 
(32) 
(33) 
(34) 



where £r = £, £2 = £3 = an d £4 = P- The pure recur- 
rence relations can be almost automatically derived. For the 
first two variables it is evidently 



£ m D m L n+1 x = £ m D m L n v = L n (f 



(35) 



(21) 



L n+1 r] = £, m D m L n+1 y = £ m D m L™w = L n p. (36) 
For the third variable one gets 
L n+1 (j> = £ m D m L n+1 v = 

— Djji, L X 

- 2 E^ UP™ [L k xL n - k y] = 
fc=0 \ / 

= -L^-2^r\ m [(D m L k x)(L n - k y) + 

k=0 V / 

+(L k x)(D m L"- k y) 

= -L^ -2^ R [(C m C m L fc x)(L- fe y) 
fe=o \ / 

+ (L k x)^ m D m L n - k y)] = 
fe=0 \ / 

+L fc x J L n - fc j ? ] . (37) 



+ 
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The same procedure can be performed for p and the result 
is 



Ln+l ttl \ I i y k j- t n — k . r k r n — k j- 

P = —L j] — I t ) |_ ^ x x 



-L k 7]L n - k y - L k yL n ' k ri I . (38) 



3 THE LIE-DERIVATIVES FOR THE iV-BODY 
PROBLEM AND ITS LINEARIZED FORM 

Let us have K point masses m 4 (i = 1, . . . , if) moving under 
the mutual gravitational attraction described by Newton's 
universal law of gravity. The coordinates and the velocities 
of these particles are denoted by x im and v im , where m is 
the index for the spatial dimension (m = 1,2,3). In the 
following sections we denote the bodies by indices i, j, k, 
. . . and the spatial indices by m, n, p, . . . , therefore the Ein- 
stein summation convention should be performed between 
the appropriate limits, which is not explicitly noted every- 
where. 

Following H&D, we present the derivation of the recur- 
rence formulae for the Lie-derivatives. The whole calcula- 
tion is presented in Appendix [A] We note that with differ- 
ent types of notations the calculation can also be found in 
H&D, some steps of the derivation should be emphasized for 
further calculations of the linearized equations. 



3.1 Equations of motion 

Using the above notations, the equations of motion of the 
JV-body problem are the following, 



Xim — ^im , 



Xim Xj 



- G ]T m] 



(39) 
(40) 



where G is the Newtonian gravitational constant and pij is 
the distance between the zth and jth body, i.e. 



Pij — ^ ^ ("^J m 3>im) — AijmAijj 



(41) 



We also introduce the following new variables and differen- 
tial operators: 





- — Xim ^jm, 


(42) 




■ — fim fjm. 


(43) 




■ — AijmBijm: 


(44) 


4>ij 


-3 

:= Pij , 


(45) 


Dim 


d 

dx im ' 


(46) 



^■in 



dVir, 



(47) 



With these notations the Lie-operator of the equations of 
motion can be written as, 



Vim D- 



G E 



/ TYlj(f>ij Aijm I Ai- 



(48) 



for the variables x im , Aij m , Bij m , Aij, v irn and 
following system of equations: 

rn+l. 



In Appendix [X] we prove that the recurrence relations 

is the 

(49) 
(50) 
(51) 

,(52) 



T n A - ■ 

-^0 Bijm 



jr n 

Ln t n 

K 



E 



j. I L (j}iiL n Ai 



J <Pij 



Pi/ FnkLg 



k cf>ijL k l Aij, 



(53) 
(54) 



where F nk = (-3)(™) + (-2)( fc ^ 1 ). We note that F nk is 
equivalent to the matrix A nk introduced in H&D. 

3.2 Linearized equations 

For the linearized coordinates and velocities we introduce 
the variables £i m and rji m , respectively. Therefore, using the 
Lemma, we get the Lie-derivatives of the linearized variables, 
namely, 



L n £im = {£,kpDk P + r]k p Ak P )L n Xi 
L n rjim = (£,k P Dkp + r]k p Akp)L rl Vi, 



(55) 
(56) 



To obtain recurrence relations we have to introduce other 
auxiliary quantities. First, we form two vectors which con- 
tain all of the linearized variables and the differential oper- 
ators: 



1 — 'kp ■ — 



£,kp 

Vkp 



(57) 
(58) 



Therefore, one can write E-V — E kp T> kp = £, kp D kp + r] k pA k p 
which simplifies the notation of the scalar products appear- 
ing in equations (I55ll - (|56l ): 



L"£, 

Ln 
1)v. 



— H ■ T>L n Vi 



(59) 
(60) 



Second, let us introduce ctij m := fi m — £ jm and f3ij m := 
Tjim — Vjm- With these newly introduced variables and ex- 
pressions we can derive the recurrence formulae for the 
linearized variables. The calculations are presented in Ap- 
pendix [B] in more details, and the result is 

L n+1 £im = L n rn, n , (61) 

L Otijm = L £im L £,jmi (62) 
L 0ijm — L 7]i m — L Tljm , (63) 

E ■ T>L n Aij — ^ I ^ J \L k a>ij m L n ~ k Bi 3m + 
fc=o \ / 

+ L k A l]m L n - k p ijm \ , (64) 



L" +1 rj im = -G ^2 m i \ E 



(65) 



Linearized Equations of the N-body Problem 5 



%i 



%5 



(S ■ VL k <t>ij)L n k Ai jm + 

i r fe j r n — k 

+L (pijL a ijrn 

n 

+p- 2 ^F nfc [% 2 ], 



(S ■ VL n ~ k <f>ij)L k kij + 
+L n - k <t> l] (E-VL k K ] ). 



(66) 



For the initialization of the recursion we have to calculate 
S • VI°<pij = 3 • V(f>ij . It is easy to show that 



■ V(/>ij = 3 • £> ft / = 



-3p,„ 5 a 



(67) 



We have some remarks concering the derivation and 
evaluation of the above formulae. First, we did not need 
the linearized equations explicitly to derive the recurrence 
relations for the linearized variables. Second, because of the 
symmetry properties of the variables, we do not have to 
calculate all of the matrix elements: we know that the ten- 
sors Aijm, Bij m , aij m , fiij m are antisymmetric for swapping 
the indices i and j and the matrices Ay, <j>ij, 3 • VAij and 
3 • V(f>ij are symmetric. Because distances are defined only 
between different bodies, the diagonal matrix elements of pa 
and their derived (cfru, L"pu, L n <f>n, 3 -T>L n pa, . . . ) are not 
defined. 



3.3 Motion in the reference frame of one of the 
bodies 

In the description of the Solar System or in perturbation 
theory, the equations of motion are transformed into a ref- 
erence frame whose origin coincides with one of the bodies. 
Practically, it is the body with the largest mass, in the So- 
lar System it is the Sun (where all orbital elements are de- 
fined relatively to the Sun) . Therefore, it could prove useful 
to have the recurrence relations both for the equations of 
motion and for the linearized part of the equations in this 
reference frame. 

Let us define the central body as the body with the 
index of i = 0. Altogether we have 1 + K bodies, where 
the other ones are indexed by i = 1, ...,K. For simplic- 
ity, denote its mass by A4 = mo. In an intertial frame, the 
equations of motion can be splitted into two parts, namely 

%im — ^im, V^^J 

(69) 
(70) 
(71) 



^im — G ^ ^ TYljpij (Xim Xjm), 



10m = V& 



in • 



K 



VOm = — G y]mjP f(xo m — Xjm)- 
2=1 

Following the usual steps, the equations of motion in the 
fixed frame can easily be derived by subtracting equa- 
tion (|7UI) from the equations of (|68|) for all i indices. Let 
us define the new variables 



Tim 
Wim 
Pi 



%im ^Om; 
pOi = Pi0, 

pT 3 - 



(72) 
(73) 



Table 1. Timing data for the Runge-Kutta integrators, the 
Bulirsch-Stoer integrator and for the Lie-integrator for different 
orders. See text for further details. 



Integrator 




otcpsize 


. [method] 

(e) toi c 




[methodj 


time 


2.4 • 10 — 


2.4 • 10 


2.4 ■ 10 


RK4 


0.302 


0.0140 


0.0091 


0.0026 


RKN5/6/ 


0.460 


0.0243 


0.0153 


0.0097 


RKN7/8/ 


0.941 


0.2521 


0.1962 


0.1548 


BS 


8.578 


2.0172 


1.7734 


1.5773 


M = 6 


0.916 


0.0829 


0.0603 


0.0437 


M = 7 


1.169 


0.1095 


0.0796 


0.0556 


M = 8 


1.421 


0.2271 


0.1776 


0.1339 


M = 9 


1.706 


0.2916 


0.2275 


0.1792 


M = 10 


2.004 


0.4332 


0.3541 


0.2912 


M = 11 


2.336 


0.5535 


0.4480 


0.3525 


M = 12 


2.689 


0.6941 


0.5762 


0.4715 


M = 13 


3.048 


0.9055 


0.7770 


0.6137 


M = 14 


3.411 


0.9352 


0.8555 


0.7113 


M = 15 


3.810 


0.9414 


0.8805 


0.8789 


M = 16 


4.223 


0.9492 


0.9414 


0.9383 



Note that the quantities pi and pij, like so cfii and 4>ij are 
distinguished only by the number of their indices. Obviously, 

Aijm — Xim Xjm — Tim Tjm and Bijm — ^im ^jm — 

Wim — Wjm- Thus using the relative (non-inertial) coordi- 
nates and velocities, the equations of motion in more com- 
pact form are 



(74) 



(75) 



Tim — Wim ) 

w^n = —G(M + mi)(j>ir 

K 

-g E m ' 

Without going into details, we preset the recurrence rela- 
tions of the Lie-derivatives, including the linearized vari- 
ables in Appendix [C] Some speed-up considerations with 
which the required number of operations can definitely be 
decreased are presented in Appendix [D] 



4 PERFORMANCE AND COMPARISONS 

We have implemented the method of Lie-integration as a 
standalone program, written in ANSI C, with the following 
capabilities. The program is able to integrate the equations 
of motion of the Af-body problem in the reference frame 
of one of the bodies (see equations (|74 p -(|75 p ) and paral- 
lelly, the program approximates the LCE by the Lyapunov 
Characteristic Indicator (LCI) of the system using the so- 
lution of the linearized equations. For the method of in- 
tegration one could use the classical fourth order Runge- 
Kutta (see lPress et al.|[l992l ) and Runge-Kutta-Ny strom in- 
tegrators (namely, RKN5/6/ and RKN7/8/, see iFehlberd 
1 19721 ; iDormand fc Princelll978l). the Bulirsch-Stoer integra- 
tor (BS, see also Press et al.lll992l ) as well the Lie-integration 
method (see equations (jC 1 p - (|C8f) and equations ()C9[) - (|C16[) 
in Appendix [Cj, up to arbitrary order M. The program is 
also able to figure out the optimal stepsizes to satisfy a pre- 
defined accuracy. The accuracy is derived using the differ- 
ences in the mean longitude which is the fastest changing 
orbital element. This type of accuracy control can be found 
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-180 -120 -60 60 120 180 -180 -120 -60 60 120 180 

Difference in mean longitude [degrees] Difference in mean longitude [degrees] 

Figure 1. The LCIs for a fictious asteroid having the same orbit as Jupiter. In the left panel, one can see the derived indicators by 
the method of RKN7/8/ (crosses) and using Lie-integration (empty squares) as the function of AA. In the right panel, the ratios x(AA) 
defined by equation 1771 1 are plotted. 



in many integrators (e.g. ORBIT9) where the dimension- 
less accuracy is defined as the difference of the mean lon- 
gitudes between the exact and approximated solution (in 
radians) divided by the the square of the number of revolu- 
tions, namely 



| AA| (radians) 
/V 2 

revolution 



(76) 



(see lMilani fc Nobilj[l988l . for a more detailed explanation). 

As an initial test, we have compared the LCIs computed 
by two different integration methods, namely RKN7/8/ and 
the Lie-integration with the order of M = 8. The dynamical 
system is the spatial Sun - Jupiter - Saturn - test parti- 
cle spatial restricted four-body system where the latter has 
the same orbit as the Jupiter has. The LCIs are calculated 
as the function of the difference in the mean longitudes of 
Jupiter (A,j) and the test particle (A m ) while all other 5 
initial orbital elements are equal to those of Jupiter. The 
results are plotted in Fig. [1] In the left panel of Fig. [1] one 
can see the derived indicators by the method of RKN7/8/, 
LCIrkn7/8/ an d using Lie-integration, LClLic as the func- 
tion of A A = A m — A j. In the right panel, the absolute 
value of the base-10 logarithm of the ratio of the indicators, 
namely 



X = 



log 



LCI 



R.KN7/8/ 



LCI L 



(77) 



are plotted resulted by this two integration method. Note 
that the integration length is 10 6 yrs, therefore the LCIs 
concerning to regular solutions are saturated around ~ 
10 -5 — 10~ 6 1/yr. It can easily be seen that in the stable re- 
gions (around the two Lagrangian points at AA = —60° and 
A A = +60°) the results of the two methods are very similar, 
the magnitude of the differences between them is « 10~ 5 . In 
the chaotic regions, the two methods yielded different LCIs 
but their magnitudes were always the same. 



tion with RK4, RKN5/6/, RKN7/8/, BS and with the Lie- 
integration and parallelly the linearized equations to get the 
result with a previously given accuracy. The ratio of the net 
CPU times is the relative cost: 



cost 



[other] 



[Lie] 
r CPU 
[other] 
r CPU 



(78) 



As one can see, the smaller the cost is the more efficient the 
Lie-integration is. It should be kept in mind that this rela- 
tive cost does not only depend on the other method but also 
on the order of the Lie- integration and the desired accuracy. 
Going into the details, the cost has been measured indirectly 
by the following way. It can be said that any of the integra- 
tion algorithms, the RK-based ones, the BS and the Lie- 
integration use the same CPU time per step independently 
from the stepsiz^E Let us denote this atomic CPU time 



by t 



[method] 
(0) 



Therefore, if the optimal stepsize At 



method] 



known for a given method and accuracy, the total CPU time 
can easily be calculated: 



[method] 
r CPU 



[method] 
r (0) 



method] : 
(<0 



At 



(79) 



where T is the total length of the integration. Because the 
relative cost is the ratio of two such value of tcpu for two 
methods, the total length of the integration cancels. The 
atomic CPU time can easily be measured, the only un- 
known is the At optimal stepsize for the different methods. 
The latter is determined by the following way. The exact 
mean longitude for the fastest rotating planet is derived for 
a given time-span (which is defined by the accuracy, see 
equation l|76[l) with an appropriately small stepsize. After 
it, the stepsize is increased iteratively by a bracketing algo- 
rithm until the integration yields a mean longitude which 
differs from the exact one by the AA value determined also 
by equation (J75J). We should note that this implies that the 



4.1 Performance analysis 

We have compared the efficiency of the Lie-integrator and 
the other implemented integrators. Here we give how much 
CPU time is required to integrate the equations of mo- 



1 In our tests, in the BS method the adaptive variation of the 
number of MMID substeps has been disabled, i.e. the extrapola- 
tion is performed after the same sequence of number of substeps: 
it yields an evaluation time which is independent from the step- 
size. 
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0.5 1 1 1 1 1 1 1 0.5 1 1 1 1 1 1 1 

7 8 9 10 11 12 13 7 8 9 10 11 12 13 

Negative logarithm of accuracy Negative logarithm of accuracy 

Figure 2. The relative cost in the CPU time as the function of the desired precision for the Lie-integration against the RKN7/8/ (left 
panel) and the Bulirsch-Stoer method (right panel) for the model system of Sun - Jupiter — Saturn. The curves show the cost for different 
orders of the Lie-integration (plus signs for M = 6, crosses: M = 8, stars: M = 10, open squares: M = 12, filled squares: M = 14 and open 
circles M = 16). The thick solid line shows the unity cost, where the RKN/BS and Lie-integration methods have the same performance. 
Smaller costs represent lower CPU usage, therefore higher gain and better performance. 




0.5 1 1 1 1 1 1 1 0.5 1 1 1 1 1 1 1 

7 8 9 10 11 12 13 7 8 9 10 11 12 13 

Negative logarithm of accuracy Negative logarithm of accuracy 

Figure 3. The relative cost in the CPU time as the function of the desired accuracy for the Lie-integration against the RKN7/8/ (left 
panel) and the Bulirsch-Stoer method (right panel) for the model system of Sun - Jupiter - Saturn - fictitious asteroid, while the linearized 
equations are also evaluated only for the massless test particle. The curves show the cost for different orders of the Lie-integration (see 
also Fig.0. 



stepsize At is constant during the integration. In Table[T]we 
summarize these timing values for some values of accuracy 
and for the Runge-Kutta methods, for the Bulirsch-Stoer 
method as well as for the Lie-integration method for orders 
M — 6, . . . , 16 for the dynamical system of Sun - Jupiter 
- Saturn - test particle extended with the linearized equa- 
tions of the latter. The second column contains the atomic 
CPU timfl while the other three columns show the opti- 
mal stepsize A ^™ cthod ', derived by the above manner for 
the accuracies e = 2.4 ■ 1CT 11 , 2.4 ■ 1CT 12 and 2.4 • 10~ 13 
respectively. Thus, the cost can be derived by the fractions 
of the appropriate values taken from this table, namely: 



[ml] A [m2] 
cost [ml]against[m2] = (0) ^je) 



2] A Jml] 



(80) 



'(()) 



"(e) 



2 measured on an Athlon XP 1800+ processor with GCCv4.1.2 
compiler, in the units of 10 — 6 seconds. 



where [ml] and [m2] index the two methods to be compared. 
We should note that timing values were not only derived 
for these values of accuracy as it can be read from Table [1] 
and we have made timing measurements when the linearized 
equations are omitted. See next sections for more details and 
for other plots. 



4.2 Efficiency as the function of the accuracy 

In Fig. [2] the relative cost of the Lie-integration against the 
RKN7/8/ and the Bulirsch-Stoer integration method are 
plotted for the three-body problem of Sun - Jupiter - Sat- 
urn as the function of the accuracy. Different curves show 
the cost for different orders of the Lie-integration between 6 
and 16. It can easily be seen that for higher orders and be- 
low a critical accuracy the Lie-integration is more efficient 
than RKN7/8/ and for higher orders, the Lie-integration 
is more efficient than the Bulirsch-Stoer method almost in- 
dependently from the accuracy. Note that in this plot the 
linearized equations are omitted from the calculations. 
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sm(x) = 



2fc+l 



(2fc + l)!' 



(81) 



10 12 
Order of the Lie-integration 

Figure 4. The relative cost in the CPU time as the function of 
the order of the Lie-integration against the RKN7/8/ method for 
the model system of Sun - Jupiter — Saturn — fictitious asteroid. 
In these runs the linearized equations are also evaluated for the 
massless test particle. The thin solid line, the long dashed line 
and the dashed line show the cost for the accuracy of 2.4 ■ 10 — 11 , 
2.4 ■ 10 -12 and 2.4 • 10 -13 , respectively. The thick line marks the 
unity cost, below which the Lie-integration is more efficient. 



In Fig. [3]we have plotted the cost of the Lie-integration 
against the methods as above but the dynamical system is 
extended with a massless test particle and for the latter 
the linearized equations are also solved. The qualitative be- 
haviour of the cost as the function of the accuracy and the 
orders of the Lie-series is almost the same as in Fig. [2] We 
note that the different methods used in the RKN, BS and 
Lie-integration to evaluate the linearized equations result 
different number of operations, therefore the costs won't be 
exactly the same. Namely, the relative CPU time cost of the 
Lie-integration against the other methods is slightly larger 
when the linearized equations are solved parallely. 

As a conclusion, we can say that omitting the linearized 
part orders below M f» 10 the Lie-integration method is infe- 
rior to the RKN7/8/, while the equations are extended with 
the linearized equations for the massless particle, the Lie- 
integration is more effective than RKN7/8/ for orders larger 
than M « 12 below a certain accuracy about e ~ 10~ . 
Comparing with the BS method, the Lie-integration is more 
effective for orders larger than M ~ 8 and M » 10 when 
the linearized equations are omitted or not, almost indepen- 
dently from the accuracy. We note that the Lie-integration 
is effective with more than a magnitude (or more) than the 
lower-order Runge-Kutta methods, as it can easily be de- 
rived from Table. [1] and equation (|80|l . 



4.3 Efficiency as the function of the order 

As it was written in the introduction, the method of Lie- 
integration approximates the Taylor-expansion of the solu- 
tion up to a finite order. One can easily prove that the appro- 
priate order, n of a Taylor-series to obtain a certain accuracy 
of a periodic function defined on an interval is proportional 
to the length, L of this interval. The concept of the proof 
is as follows. An adequately smooth periodic function can 
be approximated as a sum of sine (and cosine) functions, 
the so-called Fourier terms. The sine function, sin (a;) can be 
expanded as 



To obtain an accuracy of unity, the last (n = 2k + 1th) 
term of the series should be the solution of i° = L n m n\. 
Therefore, using Stirling's approximation, one gets 



n log L — log(n!) w n log n — n, 



(82) 



so n ~ e • L, which means n ~ L. This is true for all Fourier 
terms of the expansion of a periodic function. 

Thus one can assume that to obtain a certain accuracy, 
the M number of the terms in the Lie-series is proportional 
to the length of the integration stepsize, namely At « kM. 
The total number of arithmetical operations, therefore the 
required CPU time is a quadratic function of the order M: 
tcpu = Ordo(M 2 ). To be more precise, the required CPU 
time is tcpu = a + f3M + yM 2 , for smaller M's, the first two 
terms, a and f3M are not negligible. Therefore, to integrate 
the equations over an interval T requires 



T t , p t f} 1 «(« + /?M + 7 M 2 )-^ 



(q + 13 M + jM 2 ) 



(83) 



CPU time, which, depending on the ratios of the constants 
a, (3 and 7 has a minimum corresponding to the optimal 
order of the Lie-integration method. 

We have tested this type of dependency of the CPU time 
on the order of the Lie-integration. The results are plotted 
in Fig. [4] for three values of accuracy, while the dynamical 
system is the restricted four body problem of Sun - Jupiter 
- Saturn - test particle extended with the linearized equa- 
tions respecting to the latter. As it was assumed above, the 
relative cost has a minimum corresponding to the optimal 
order and the value of this minimum is around M ~ 13 — 15, 
depending on the accuracy. It can also be seen that for larger 
values of M the cost increases which means worse efficiency, 
as it is expected from equation (|83p . What is more interest- 
ing that the position of the minimum clearly depends on the 
accuracy: the better the accuracy is the larger the optimal 
value of M is. This might imply another kind of adaptive 
integration method where not only the stepsize varies but 
the order of the Lie-integration. 



4.4 Implemetaion of the method 

We have implemented the method of Lie-integration of the 
JV-body problem as a standalone ANSI C program, extended 
with the linearized equations and the capability to calculate 
the LCIs. The version of the program which can integrate the 
motion of 1+3 bodies and was used in our benchmarks can 
be downloaded from the address http : //cm . elte . hu/lie as 
a single .tar.gz archive. The full version which is capable 
to integrate the motion of arbitrary number of bodies can 
be requested from the first author via e-mail. All versions of 
this code are designed to work on UNIX-like environments. 



5 DISCUSSION AND SUMMARY 

In this paper we have presented a lemma with which re- 
currence relations can be derived for the Lie-integration of 
linearized equations. We have demonstrated the usage of 
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this lemma on the Henon-Heiles system, and thereafter ap- 
plied it to the equations of the iV-body problem, includ- 
ing the non-inertial equations where the origin of the refer- 
ence frame is fixed to one of the bodies. Our performance 
comparisons have shown that although these recurrence for- 
mulae are rather complicated, they can efficiently be used 
for integrations where high accuracy is required. We have 
investigated realistic dynamical systems for these compar- 
isons: using the lemma, the recurrence relations were deter- 
mined and using the Lie-integration technique, the LCIs for 
a fictitious asteroid was computed. The method of LCIs is 
the basis for many modern chaos detection methods, there- 
fore our lemma and the derived Lie-integration method can 
widely be used in various kind of dynamical investigations, 
providing a faster alternative to the currently used tech- 
niques. We have checked the efficiency as the function of 
the order and accuracy. These tests have shown that the 
Lie-integration is definitely more effective than the classical 
RK4 and RKN5/6/ integration method and above a certain 
accuracy about e ~ ICP , the Lie-integration is more ef- 
fective than the method of RKN7/8/ for orders larger than 
M — 1 or M — 8, whether the linearized equations are eval- 
uated parallelly or not. We found that the Lie-integration is 
more effective than the BS integrator for orders larger than 
M ~ 10, almost independently from the accuracy. 

Further studies are already ongoing concerning this 
problem. First, there could be several possibilities for opti- 
mization in the actual implementation of the Lie- integration: 
we expect that the re-ordering of the highly nested loops 
and/or the introduction of new auxiliary variables yield bet- 
ter performance. Second, some aspects of the Lie-integration 
should better be analysed and understood, including the 
long-term error propagation in higher orders which is the 
basis of the adaptive extensions of this integration method. 
And last, we are going to develop a more general code which 
is not only capable of the calculation of LCIs but can be ex- 
tended with other chaos indicators. 
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APPENDIX A: RECURRENCE RELATIONS FOR THE TV-BODY PROBLEM 



Using the notations denned in Section 13.11 we derive the recurrence relations for the equations of motion for the TV-body 
problem. 

As we have shown in Section ^. li the Lie-operator of the equations of motion is 



Lq — vi m Di m — g y ' 



(Al) 



This implies that the first Lie-derivatives of the coordinates and velocities are the right-hand side of the equations of motion, 
namely 



LoXim — v., 

LoVim — — 



(A2) 



G ^ ^ UljPij {Xjm ^im) — G ^ TYlj(j)ij Aij m . 



The distance py does not depend on the velocities, like so (j>ij, therefore their Lie-derivatives can easily be calculated: 



L(3pij = VkmD, 



kmPij • 



(A3) 



Therefore, 



Lopij 



E 



E 



l-'kr. 



dx k 



I'kn 



2pij 



(2(Xjm — Xim)(5kj ~ Ski)) 



Y m 

'^2 v kmPTj 1 [(Xjm ~ X im )(Skj — Ski)] = 



~ Pij [("jm H, m )(lj m Xi m )] — p^ Bji m Aji m — /5y Aji — Ai j . 

Now one can calculate the Lie-derivative of cpij = Pij 3 '- 

L (f>ij = Loip'/) = -ip~ A L pij = -3p~/Aij = p~ 2 (-?><j)ijAij). 

With mathematical induction, one can prove that 



m+l 



'V = Pi/ E 



- 2 



II 

k + 1 



j i^ij-^O^ij • 



(A4) 
(A5) 

(A6) 



For n = this equation is equivalent to equation (|A5[) . Let us assume that this is true for all m ^ n, and calculate Lq +2 0,j: 



b H = L o ft/ E 



- 2 



ii 

k + 1 



= L (p i:j 2 ) ■ (p 2 jLo +1 (jjij) + pi/ ^2 



- 2 



n 

k + 1 



Lq^ (j)ij Lq Aij + 



+ft/E 



k=0 L 



k + 1 



Lq k (j)ij Lq + 1 Aij — %. 



(A7) 



The first term is 

L ° (Pi/) ' {Pij L o +1 ^n) = - 2 Pij 1 ^j ( PijPij L o +1 <l ) ij = Pi? (-2io + VijA i: ,) 



(A8) 



We can increase the upper limit of the first summation of term 2 in equation (|A7j) from n to n+1, since in the appearing new 
terms, the factors ( n ? 1 ) an 
in term 3 of equation l|A7fl : 



terms, the factors ( n ? 1 ) and ( n " 2 ) are zero by definition. To unify term 1 and term 3, we introduce a new index, k' = k + 1 



: E 



k + 1 



n+l 

T n-k i T k+1 . -2 ST^ 

L <pijL Aij = p^ 2^ 



k' - 1 



- 2 



J^Q tyii-L'O H-ij. 



(A9) 



Note that if we substitute k! = into the expression after the summation, we get the same what equation (|A8[) is, therefore 
the latter can be inserted into the summation of equation (|A9f) while the lower limit of k' = 1 is replaced to k! — 0. Therefore: 
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n+l 

% = pit E 

fc'=0 



Using the relation (£) + ( fe ^) = we get 



^HU^HU-J-^MT Hill 



therefore we could simplify equation 1A10[) : 



n + l 
fc=0 



t n+l — k i j k * 

-^0 (pij-l^o-lH' 



Comparing equation (|A12j) with equation (| A6|) . we conclude that the relation is proven. For simplicity, we define 



k + 1 



(A10) 



(All) 



(A12) 



(A13) 



Continuing the derivation of the recurrence formulae, we calculate the higher order Lie-derivatives of Aij using the binomial 
theorem: 

n / \ n / \ 

(A14) 



Lq A.ij — ^ | - J Z/Q-Aij m Z/Q Bijm — S ^ I J (LtQXim ~ ^0 %jm) {Lq Vim Lq ^. 



In the equations of motion the term cf>ijAijm appears, its higher order Lie-derivatives can also be calculated like the last 
relation for LqA^: 



Lq {4>ijAij m ) — S ^ I jZ/o^ij {Lq %im — Lq Xj m ). 



(A15) 



n I rfc 

To summarize our results the complete set of the recurrence relations for the equations of motion can be found in equations ((49 
(|54l). 



APPENDIX B: DERIVATION OF THE LINEARIZED EQUATIONS 

To obtain the recurrence relations for the linearized equations, we apply the operator S • 23 to equations (|49I) - H54I ). We note 
that the operator E • 23 is linear, 



S ■ V(pa + qb) = p(E ■ 23a) + q(3 ■ 23b), 

where a and 6 are continuous functions while p and q are constants, and one can use Leibniz's rule: 
S ■ V(ab) = (S ■ Va)b + o(H ■ Vb). 

Moreover, we should note that the operators S ■ 23 and L cannot be commuted, E • 23 L 7^ LS ■ 23. 
For the first three equations, we get 

Ln+1> 1 — 1 -T-iT-n + 1 1 — 1 -pi r n _ j n 

C,im — " ' ^im — >— ' " Ul-J Vim — i-t 7)im, 

Ln r n /- T n C 

OLijm — i-i Sim J-J sj'm, 

Pijm ^ T]zm Li TJj m . 

For S • 23 Aij we can use the linear property and apply Leibniz's rule: 



■ VL n A t] = 



f . J |^(S • 23L fc Aij m )L n k Bij m + L k Aij m (E ■ VL n fc _B i:)m )j 
h=o \ ) 

n / \ 

^ ^ I J ^3 rn ^ J Bij m ~t~ L AijmL /?ijm)J . 



(Bl) 
(B2) 



(B3) 
(B4) 
(B5) 



(B6) 



Here we applied the identities E ■ T>L n Aij m 
follow the same procedure: 



L n ctijm and S ■ T>L n Bij m = L n f3ij m - For the calculation of E • T>L n+1 (f>ij 



S-23L" 



(S-23) 



' F n hL n k <pijL k Aij 



n 

(E ■ Vp^)(p%L n+1 ^) + E F nk [(S ■ VL n -*4> ii )L k A ii + ^"^(S ■ 23L fe Ay)] 



(B7) 
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The only unknown factor in equation (|B7|I is the quantity S ■ T>p^ . We can calculate it easily, because p t j only depends by 
definition on the coordinates: 

S ■ £>Pi/ — £,kmDkm(pij 2 ) — ^km{~2)p i:j 3 DkmPij = ( — 2p^ 3 ) £km.DkmPij = %• (B8) 

The expression ^kmDkmpij can be calculated like equation (|A3[l . where we replace Vkm by £fc m : 

^km-^kmPij ~ Pij (^im ^jm)(^im = Pjj O^ijm 

A ljm . (B9) 

Thus, we get 

% = -2py 4 ay m Ay m . (BIO) 
Adding all terms together, we obtain 

n 

H ■ = -2pr?a ijm A ijm L n+1 <t> i:j + P^Y, Fnk [( H ' VL n ~ k <kj)L k Mi + i n ~ fc ^«(S ■ M, fc Ay)] . (Bll) 

The derivation of L n+1 r/i m = S • VL n+1 Vim is the following: 

L" +1 r/ lm = S • VL n+1 v im — — G ^ m J ^ M 5 ■ D ^L k </> ij L n - k A ijm ) \ = 

3=l,J^i U=0 V / J 

= -C ^ mJf^ ("J [(H • M^y)!"-^ + (S • M n - fe 4 ijra )] 
i=i,j& U=o V / 

= -G J] m ife U) [(S-^^yJi^-Atim+i^i^oym]}. (B12) 

3 = 1, J^i (k=0 V / J 

Now we obtained the recurrence relations for all of the linearized coordinates, velocities and the auxiliary variables 5 • 2? Ay , 



S • T>(f>ij. The complete set of these equations are summarized in equations (|61[) - (|66[ ). 

APPENDIX C: MOTION IN A REFERENCE FRAME FIXED TO ONE OF THE BODIES 

Throughout the derivation of the recurrence relations, we can use the fact that the partial differential operators g ^ and — 

are equivalent to Di m — — and Ai m = — , because the variables differ only in a constant (xo m and Vo m , respectively). 
We have to define the new variable Aj = ri m Wi m . The derivation of the recurrence relations can be done following the steps 
of Appendix [X] and Appendix [B] the quantities p;, cpi and Aj have the same properties for the Lie-derivation as pij, <j)ij and 
Ay , respectively, therefore all of the induction steps can be done in the appropriate way. 

Thus, the recurrence relations for the Af-body problem around a fixed centre can be written as 

L n+l r im = L n w im , (CI) 

L Aij m = L Tim L Tjm, (C2) 
L Bijrn = L Wi m — L Wjrn, (C3) 

L n A t = J2\l) Lkr ^ L "~ kw ^ (° 4 ) 



fc=0 



-G{M + mi ) U J L k faL n - k r im -G m ^J2 (I J [L k 4>^L n ' k A ijm + L k <j> J L n - k r jm ] , (C6) 

fe = V / j=l,j^i k=0 \ / 

n 

! J2PnkL n - k 4> i L k A i , (C7) 

fc=0 
n 

L n+1 &3 = p- 2 ^F nfe L n - fc ^L fc Ay. (C8) 

k=0 

Let us denote the linearized of r im and W; m by £i m and rj irn , respectively. Since Uij m = £i m — £j m and /3y m = r\ im — rjj m , for 
the linearized equations the calculations yield 

L n+1 £ lm = L n r) im , (C9) 

L OLijm L (im L (j m , (CIO) 

L /3ijm L Tjim L Tjjm, (Cll) 
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3 ■ VL n Ai = ^2 f ^ ) { Lk ^imL n ~ k w im + L k r im L n ~ k i] lm }j , (C12) 

fc=0 V / 

E-T>L n Aij = 1 J {h h ctij m L n k Bij m + L k Aij m L n fc /3 i:)m ^ , (013) 

L n+1 r, im = -G(«M + m.) U J [( H ■ VL k <f>i)L n - k r im + L k ^L n ~ k ^ - 

K n / \ 

-G ^ m J ^ I ? J [(- ■ r DL k <t> ij )L n ~ k A i]m + L k (j>ijL n ~ k ai jrn + (3 ■ VL k (t> ] )L n ~ k r ]m + L k (f>jL"~ k $fii^) 
: ■ M n+1 & = -2pr 2 ^mr lm L" +1 ^ + ft 2 I] f»* [(3 ■ T>L n - h fc)L k As + L"- fc ^(H ■ OL^A,)] , (015) 

n 

■ VL n+1 (f)ij = -2pr 2 Q, jmy 4 ljm L n+1 ^ +Py 2 5Z^ [(S-ZJL"-*0y)L*A y +L"-*^ y (S • 2?L fc A y )] . (016) 



APPENDIX D: SPEED-UP CONSIDERATIONS 

Introducing new variables, the required number of arithmetical operations can be decreased in equations (|C1[) - (1C8[1 and 
equations (|C9|I - (|C16|) . Namely, the calculation of L Wi m and L Tjim can be written as 

L" +1 w m = -G(M + mi )sW-G J2 mi(sW m + sW), (Dl) 

L n+1 Vtm = -GiM + m^W-G jh m,(EW +EW), (D2) 
where the new variables are 



£ = 

zm 

fc=0 



S ™ = E (D3) 

fc=0 V / 
fe=0 V / 

jh ( ?) [(3 ■ VL k d>i)L n - k r im + , (D5) 

E ( j! ) [( S ' 'DL k (f>ij)L n ~ k A i: j m + L k (t>i ] L n ~ k o> l j rr ^ . (D6) 

The implementation of the above relations can increase the speed of the calculations by 20%-30%, depending on the number 
of the bodies. 
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